Evolvability

In the present report, we first estimate the narrow-sense heritability and evolvability at the species level, i.e. across all populations. Then we estimate these two parameters at the population level, i.e. one parameter value per population. For that, we first run one model per population (hereafter referred as ‘population-specific models’). Second, we run a full model with all populations (hereafter referred as ‘full model’), based on Archambeau et al. (2023).

As heritability and evolvability are specific to the environment, the models are fitted in each common garden separately. We fit models that do not account for the nursery effect for the two common gardens in which trees come from the same nursery: the sites of Yair (FS) in the Scottish borders (in which all trees come from NG) and Glensaugh (FE) in which all trees come from NE (except four trees that come from NG and that we remove for the analyses).

We fit models that do account for the nursery effect in the common garden in which trees come from different nurseries: Inverewe (FW), which contains cohorts of trees raised in each of the three nurseries as follows: 290 grown locally in the NW; 132 grown in the NG; and 82 grown in the NE.

The data

For details on the data, see report EstimatingPhenotypicPlasticityScotsPine.qmd.

Code

# We load the field data
data <- read_excel(here("data/Field.xlsx"), na = "NA")

# we load the nursery data (dataset updated by Annika - 02/06/2023)
nursery_data <- read_excel(here("data/SPnurseries_updatedByAnnika02062023.xlsx")) %>% 
  dplyr::rename(FieldCode=Tag)

# we merge the two
data <- data %>% left_join(nursery_data, by=c("FieldCode")) 

We have to change the names of the block to specify that they are nested within sites:

Code

# site-specific blocks
data <- data %>% mutate(Block = paste0(FieldSite,"_",Block))

According to Beaton et al. (2022), there are four randomised blocks in both FS and FE and three in FW. We look at the number of individuals in each block:

Code

data %>% group_by(Block) %>% summarise("Number of individuals"=n()) %>% kable_mydf()

Block	Number of individuals
FE_A	168
FE_B	168
FE_C	168
FE_D	168
FS_A	168
FS_B	168
FS_C	168
FS_D	168
FW_A	168
FW_B	168
FW_c	1
FW_C	167

We see that there is one typo in FW (block noted as c instead of C), that we have to correct:

Code

data <- data %>% mutate(Block=case_when(Block=="FW_c" ~ "FW_C",
                                TRUE ~ Block))

data %>% group_by(Block) %>% summarise("Number of individuals"=n()) %>% kable_mydf()

Block	Number of individuals
FE_A	168
FE_B	168
FE_C	168
FE_D	168
FS_A	168
FS_B	168
FS_C	168
FS_D	168
FW_A	168
FW_B	168
FW_C	168

In the present report, the estimation of heritability and evolvability is done for:

Code

trait <- "HA20"

We subset the dataset by removing missing values for the trait considered and by removing the four trees in FE that belong to the NG nursery:

Code

site_codes <- unique(data$FieldSite)

subdata <- data %>% 
  dplyr::select(FieldSite,Block,PopulationCode,Family,Nursery,all_of(trait)) %>% 
  drop_na() %>% 
  dplyr::rename(site=FieldSite,
                bloc=Block,
                pop=PopulationCode,
                fam=Family,
                nurs=Nursery,
                y=any_of(trait)) %>% 
  filter(!(site=="FE" & nurs =="NG")) # we remove the four trees that come from the nursery NG in the field site FE

Species-level parameters

In this section, we estimate evolvability (\(I\)) and narrow-sense heritability (\(h^2\)) at the species level, i.e. one value for all populations.

Model equation

Here is the mathematical model of the model for a common garden with different nurseries. The mathematical model of the model for a common garden with trees coming from one nursery is the same, without the nursery intercepts.

We model each trait \(y\) such as:

\[\begin{equation} \begin{aligned} y_{bpfnr} & \sim \mathcal{N}(\mu_{bpfn},\sigma^{2}_{r})\\[3pt] \mu_{bpc} & = \beta_{0} + B_{b} + P_{p} + F_{f} + N_{n} \\[3pt] \beta_{0} & \sim \mathcal{N}(\mu_{y},2)\\[3pt] \begin{bmatrix} B_{b} \\ P_{p} \\ F_{f} \\ N_{n} \end{bmatrix} & \sim \mathcal{N}\left(0, \begin{bmatrix} \sigma^{2}_{B}\\[3pt] \sigma^{2}_{P}\\[3pt] \sigma^{2}_{F}\\[3pt] \sigma^{2}_{N}\\[3pt] \end{bmatrix} \right)\\ \end{aligned} \end{equation}\]

\(\beta_{0}\) is the global intercept. \(B_{b}\) are the block intercepts with a variance \(\sigma^{2}_{B}\). \(P_{p}\) are the population intercepts with a variance \(\sigma^{2}_{P}\). \(F_{f}\) are the family intercepts with a variance \(\sigma^{2}_{F}\). \(N_{n}\) are the nursery intercepts with a variance \(\sigma^{2}_{N}\). \(\sigma^{2}_{r}\) is the residual variance.

Partitioning of the total variance \(\sigma_{tot}^{2}\):

\[\begin{equation} \begin{aligned} \sigma_{tot}^{2} & = \sigma_{r}^{2} + \sigma_{B}^{2} + \sigma_{P}^{2} + \sigma_{F}^{2} + \sigma_{N}^{2} \\[3pt] \sigma_{r} & = \sigma_{tot} \times \sqrt(\pi_{r})\\[3pt] \sigma_{B} & = \sigma_{tot} \times \sqrt(\pi_{B})\\[3pt] \sigma_{P} & = \sigma_{tot} \times \sqrt(\pi_{P})\\[3pt] \sigma_{F} & = \sigma_{tot} \times \sqrt(\pi_{F})\\[3pt] \sigma_{N} & = \sigma_{tot} \times \sqrt(\pi_{N})\\[3pt] \sigma_{tot} & \sim \mathcal{S}^{*}(0,1,3) \end{aligned} \end{equation}\]

with \(\sum_{l}^{5}\pi_{l}=1\).

Here is the Stan code of the model:

Code

model_different_nurseries <- stan_model(here("scripts/stan_models/StandardQuantGenModel_3.stan"))
model_one_nursery <- stan_model(here("scripts/stan_models/StandardQuantGenModel_4.stan"))
model_different_nurseries

S4 class stanmodel 'anon_model' coded as follows:
// Variance partitioning in the common gardens with trees from different nurseries
data {
  int<lower=1> n;                                                               // Number of observations
  vector[n] y;                                                                  // Response variable
  int<lower=0> n_bloc;                                                          // Number of blocks
  int<lower=0> n_pop;                                                           // Number of populations
  int<lower=0> n_fam;                                                           // Number of families
  int<lower=0> n_nurs;                                                          // Number of nurseries
  int<lower=0, upper=n_bloc> bloc[n];                                           // Blocks
  int<lower=0, upper=n_pop> pop[n];                                             // Populations
  int<lower=0, upper=n_fam> fam[n];                                             // Families
  int<lower=0, upper=n_nurs> nurs[n];                                           // Nurseries
}
parameters {
  real beta0; // global intercept
  simplex[5] pi;
  real<lower = 0> sigma_tot;

  vector[n_bloc] z_block;
  vector[n_pop] z_pop;
  vector[n_fam] z_fam;
  vector[n_nurs] z_nurs;

}
transformed parameters {
  real R_squared;

  real<lower = 0>  sigma_r;
  real<lower = 0>  sigma_block;
  real<lower = 0>  sigma_pop;
  real<lower = 0>  sigma_fam;
  real<lower = 0>  sigma_nurs;

  vector[n_bloc] alpha_block;
  vector[n_pop] alpha_pop;
  vector[n_fam] alpha_fam;
  vector[n_nurs] alpha_nurs;

  vector[n] mu; // linear predictor

  // variance partitioning with the simplex pi
  sigma_r = sqrt(pi[1]) * sigma_tot;
  sigma_block = sqrt(pi[2]) * sigma_tot;
  sigma_pop = sqrt(pi[3]) * sigma_tot;
  sigma_fam= sqrt(pi[4]) * sigma_tot;
  sigma_nurs= sqrt(pi[5]) * sigma_tot;

  alpha_block = z_block*sigma_block;
  alpha_pop = z_pop*sigma_pop;
  alpha_fam = z_fam*sigma_fam;
  alpha_nurs = z_nurs*sigma_nurs;

  mu = rep_vector(beta0, n) + alpha_block[bloc] + alpha_pop[pop] + alpha_fam[fam] + alpha_nurs[nurs];
  R_squared = 1 - variance(y - mu) / variance(y);
}
model{
  //Priors
  beta0 ~ normal(mean(y),20);
  sigma_tot ~ student_t(3, 0.0, 1.0);

  z_block ~ std_normal();
  z_pop ~ std_normal();
  z_fam ~ std_normal();
  z_nurs ~ std_normal();

  // Likelihood
  y ~ normal(mu, sigma_r);
}
generated quantities {
  //Variances, narrow-sense heritability and evolvability
  real<lower=0> sigma2_tot;
  real<lower=0> sigma2_r;
  real<lower=0> sigma2_block;
  real<lower=0> sigma2_pop;
  real<lower=0> sigma2_fam;
  real<lower=0> sigma2_nurs;
  real<lower=0> h2;
  real<lower=0> I;

  sigma2_tot = square(sigma_tot);
  sigma2_r = square(sigma_r);
  sigma2_block = square(sigma_block);
  sigma2_pop = square(sigma_pop);
  sigma2_fam = square(sigma_fam);
  sigma2_nurs = square(sigma_nurs);
  h2 = sigma2_fam/(sigma2_r+sigma2_fam);
  I = sigma2_fam/(mean(y)^2);

} 

Running the models

Options for the Stan models:

Code

# Sampling in Bayesian models
n_chains <- 4 # number of chains (MCMC)
n_iter <- 2500 # number of iterations
n_warm <- 1250 # number of iterations in the warm-up phase
n_thin <- 1 # thinning interval
save_warmup <- FALSE 

# Credible intervals
conf_level <- 0.95

Code

lapply(site_codes, function(site_code){
  
  subsite <- subdata %>% 
    filter(site==site_code)
  
  if(site_code %in% c("FE","FS")){ 
    model <- model_one_nursery
    pars <- c("beta0", "pi", "R_squared", "h2","I",
              "alpha_block","alpha_pop","alpha_fam",
              "sigma2_r","sigma2_block","sigma2_pop","sigma2_fam","sigma2_tot")
  } else if(site_code=="FW"){
    model <- model_different_nurseries
    pars <- c("beta0", "pi", "R_squared", "h2","I",
              "alpha_block","alpha_pop","alpha_fam","alpha_nurs",
              "sigma2_r","sigma2_block","sigma2_pop","sigma2_fam","sigma2_nurs","sigma2_tot")
      
    }
  
  stanfit <- sampling(model,
                    data = compose_data(subsite),
                    pars=pars,
                    iter = n_iter, 
                    chains = n_chains, 
                    cores = n_chains,
                    save_warmup = save_warmup,
                    thin=n_thin)
}) %>% 
  setNames(site_codes) %>% 
  saveRDS(file=here(paste0("outputs/Evolvability/",trait,"_SpeciesLevelModels.rds")))

Warning: There were 3 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Interval plots

Variance partitioning

Code

readRDS(file=here(paste0("outputs/Evolvability/",trait,"_SpeciesLevelModels.rds"))) %>% 
  lapply(function(x){
    
  broom.mixed::tidyMCMC(x,pars=("pi"),
                droppars = NULL, 
                ess = F, rhat = F, 
                conf.int = T, conf.level = conf_level)  
  }) %>% 
  bind_rows(.id = "site") %>% 
  mutate(prop_var = case_when(term == "pi[1]" ~ "Residuals",
                              term == "pi[2]" ~ "Blocks",
                              term == "pi[3]" ~ "Populations",
                              term == "pi[4]" ~ "Families",
                              term == "pi[5]" ~ "Nurseries")) %>% 
  ggplot(aes(y = prop_var, x = estimate, xmin = conf.low, xmax = conf.high, color=site)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  ylab("") +
  xlab("Proportion of variance explained") +
  labs(color = "Field sites") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.y=element_blank(),
        panel.grid.major.y=element_blank())  +
    guides(color=guide_legend(ncol=1))

Narrow-sense heritability

Code

readRDS(file=here(paste0("outputs/Evolvability/",trait,"_SpeciesLevelModels.rds"))) %>% 
  lapply(function(x){
    
  broom.mixed::tidyMCMC(x,pars=("h2"),
                droppars = NULL, 
                ess = F, rhat = F, 
                conf.int = T,conf.level = conf_level)  
  }) %>% 
  bind_rows(.id = "site") %>% 
  
  ggplot(aes(y = site, x = estimate, xmin = conf.low, xmax = conf.high)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  ylab("") +
  xlab("Narrow-sense heritability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.y=element_blank(),
        panel.grid.major.y=element_blank())  +
    guides(color=guide_legend(ncol=1))

Evolvability

Code

readRDS(file=here(paste0("outputs/Evolvability/",trait,"_SpeciesLevelModels.rds"))) %>% 
  lapply(function(x){
    
  broom.mixed::tidyMCMC(x,pars=("I"),
                droppars = NULL, 
                ess = F, rhat = F, 
                conf.int = T,conf.level = conf_level)  
  }) %>% 
  bind_rows(.id = "site") %>% 
  
  ggplot(aes(y = site, x = estimate, xmin = conf.low, xmax = conf.high)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  ylab("") +
  xlab("Evolvability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.y=element_blank(),
        panel.grid.major.y=element_blank())  +
    guides(color=guide_legend(ncol=1))

Population-specific parameters

In this section, we estimateevolvability (\(I\)) and narrow-sense heritability (\(h^2\)) at the population level, i.e. one value for each population.

Full dataset

We first fit the models on the entire dataset.

Population-specific models

We first fit one model per population.

Here is the model equation for a common garden with different nurseries. The model equation for a common garden with trees coming from one nursery is the same, without the nursery intercepts.

We model each trait \(y\) such as:

\[\begin{equation} \begin{aligned} y_{bfnr} & \sim \mathcal{N}(\mu_{bfn},\sigma^{2}_{r})\\[3pt] \mu_{bpc} & = \beta_{0} + B_{b} + F_{f} + N_{n} \\[3pt] \beta_{0} & \sim \mathcal{N}(\mu_{y},2)\\[3pt] \begin{bmatrix} B_{b} \\ F_{f} \\ N_{n} \end{bmatrix} & \sim \mathcal{N}\left(0, \begin{bmatrix} \sigma^{2}_{B}\\[3pt] \sigma^{2}_{F}\\[3pt] \sigma^{2}_{N}\\[3pt] \end{bmatrix} \right)\\ \end{aligned} \end{equation}\]

\(\beta_{0}\) is the global intercept. \(B_{b}\) are the block intercepts with a variance \(\sigma^{2}_{B}\). \(F_{f}\) are the family intercepts with a variance \(\sigma^{2}_{F}\). \(N_{n}\) are the nursery intercepts with a variance \(\sigma^{2}_{N}\). \(\sigma^{2}_{r}\) is the residual variance.

Partitioning of the total variance \(\sigma_{tot}^{2}\):

\[\begin{equation} \begin{aligned} \sigma_{tot}^{2} & = \sigma_{r}^{2} + \sigma_{B}^{2} + \sigma_{F}^{2} + \sigma_{N}^{2} \\[3pt] \sigma_{r} & = \sigma_{tot} \times \sqrt(\pi_{r})\\[3pt] \sigma_{B} & = \sigma_{tot} \times \sqrt(\pi_{B})\\[3pt] \sigma_{F} & = \sigma_{tot} \times \sqrt(\pi_{F})\\[3pt] \sigma_{N} & = \sigma_{tot} \times \sqrt(\pi_{N})\\[3pt] \sigma_{tot} & \sim \mathcal{S}^{*}(0,1,3) \end{aligned} \end{equation}\]

with \(\sum_{l}^{4}\pi_{l}=1\).

Stan code:

Code

model_pop_specific_different_nurseries <- stan_model(here("scripts/stan_models/StandardQuantGenModel_5.stan"))
model_pop_specific_one_nursery <- stan_model(here("scripts/stan_models/StandardQuantGenModel_6.stan"))
model_pop_specific_different_nurseries

S4 class stanmodel 'anon_model' coded as follows:
// Variance partitioning in one population in the common gardens with trees from different nurseries
data {
  int<lower=1> n;                                                               // Number of observations
  vector[n] y;                                                                  // Response variable
  int<lower=0> n_bloc;                                                          // Number of blocks
  int<lower=0> n_fam;                                                           // Number of families
  int<lower=0> n_nurs;                                                          // Number of nurseries
  int<lower=0, upper=n_bloc> bloc[n];                                           // Blocks
  int<lower=0, upper=n_fam> fam[n];                                             // Families
  int<lower=0, upper=n_nurs> nurs[n];                                           // Nurseries
}
parameters {
  real beta0; // global intercept
  simplex[4] pi;
  real<lower = 0> sigma_tot;

  vector[n_bloc] z_block;
  vector[n_fam] z_fam;
  vector[n_nurs] z_nurs;

}
transformed parameters {
  real R_squared;

  real<lower = 0>  sigma_r;
  real<lower = 0>  sigma_block;
  real<lower = 0>  sigma_fam;
  real<lower = 0>  sigma_nurs;

  vector[n_bloc] alpha_block;
  vector[n_fam] alpha_fam;
  vector[n_nurs] alpha_nurs;

  vector[n] mu; // linear predictor

  // variance partitioning with the simplex pi
  sigma_r = sqrt(pi[1]) * sigma_tot;
  sigma_block = sqrt(pi[2]) * sigma_tot;
  sigma_fam= sqrt(pi[3]) * sigma_tot;
  sigma_nurs= sqrt(pi[4]) * sigma_tot;

  alpha_block = z_block*sigma_block;
  alpha_fam = z_fam*sigma_fam;
  alpha_nurs = z_nurs*sigma_nurs;

  mu = rep_vector(beta0, n) + alpha_block[bloc] + alpha_fam[fam] + alpha_nurs[nurs];
  R_squared = 1 - variance(y - mu) / variance(y);
}
model{
  //Priors
  beta0 ~ normal(mean(y),20);
  sigma_tot ~ student_t(3, 0.0, 1.0);

  z_block ~ std_normal();
  z_fam ~ std_normal();
  z_nurs ~ std_normal();

  // Likelihood
  y ~ normal(mu, sigma_r);
}
generated quantities {
  //Variances, narrow-sense heritability and evolvability
  real<lower=0> sigma2_tot;
  real<lower=0> sigma2_r;
  real<lower=0> sigma2_block;
  real<lower=0> sigma2_fam;
  real<lower=0> sigma2_nurs;
  real<lower=0> h2;
  real<lower=0> I;

  sigma2_tot = square(sigma_tot);
  sigma2_r = square(sigma_r);
  sigma2_block = square(sigma_block);
  sigma2_fam = square(sigma_fam);
  sigma2_nurs = square(sigma_nurs);
  h2 = sigma2_fam/(sigma2_r+sigma2_fam);
  I = sigma2_fam/(mean(y)^2);

} 

We run the models:

Code

lapply(site_codes, function(site_code){
  
 subsite <- subdata %>% 
    filter(site==site_code) %>% 
    arrange(pop)
  
 lapply(unique(subsite$pop), function(pop_i){
   
   subpop <- subsite %>% 
    filter(pop==pop_i) 
   
   if(site_code %in% c("FE","FS")){ 
    model <- model_pop_specific_one_nursery
    pars <- c("beta0", "pi", "R_squared", "h2","I",
              "alpha_block","alpha_fam",
              "sigma2_r","sigma2_block","sigma2_fam","sigma2_tot")
  } else if(site_code=="FW"){
    model <- model_pop_specific_different_nurseries
    pars <- c("beta0", "pi", "R_squared", "h2","I",
              "alpha_block","alpha_fam","alpha_nurs",
              "sigma2_r","sigma2_block","sigma2_fam","sigma2_nurs","sigma2_tot")
       }

 
  stanfit <- sampling(model,
                      data = compose_data(subpop),
                      pars=pars,
                      iter = n_iter, 
                      chains = n_chains, 
                      cores = n_chains,
                      save_warmup = save_warmup,
                      thin=n_thin)   
 }) %>% 
   setNames(unique(subsite$pop))
}) %>% 
  setNames(site_codes) %>% 
  saveRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels.rds")))

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 10 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 7 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 2 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 3 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Interval plots

Narrow-sense heritability

Code

df <- readRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels.rds"))) %>% 
  lapply(function(site_i){
    
    site_i %>% lapply(function(pop_i){
    
  broom.mixed::tidyMCMC(pop_i,pars=("h2"),
                droppars = NULL, 
                ess = F, rhat = F, 
                conf.int = T,conf.level = conf_level) 
      
    }) %>% bind_rows(.id = "pop")
  }) %>% bind_rows(.id = "site")


  
# df %>%   
#   ggplot(aes(y = pop, x = estimate, xmin = conf.low, xmax = conf.high,color=site)) +
#   geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
#   ylab("") +
#   xlab("Narrow-sense heritability") +
#   theme(axis.text = element_text(size=18),
#         panel.grid.minor.y=element_blank(),
#         panel.grid.major.y=element_blank())  +
#     guides(color=guide_legend(ncol=1))


df %>%   
  ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  xlab("") +
  ylab("Narrow-sense heritability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank())  +
    guides(color=guide_legend(ncol=1))

Evolvability

Code

df <- readRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels.rds"))) %>% 
  lapply(function(site_i){
    
    site_i %>% lapply(function(pop_i){
    
  broom.mixed::tidyMCMC(pop_i,pars=("I"),
                droppars = NULL, 
                ess = F, rhat = F, 
                conf.int = T,conf.level = conf_level) 
      
    }) %>% bind_rows(.id = "pop")
  }) %>% bind_rows(.id = "site")
  
# 
# df %>% ggplot(aes(y = pop, x = estimate, xmin = conf.low, xmax = conf.high,color=site)) +
#   geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
#   ylab("") +
#   xlab("Evolvability") +
#   theme(axis.text = element_text(size=18),
#         panel.grid.minor.y=element_blank(),
#         panel.grid.major.y=element_blank())  +
#     guides(color=guide_legend(ncol=1))


df %>%  ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  xlab("") +
  ylab("Evolvability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank())  +
    guides(color=guide_legend(ncol=1))

Full model

Here is the model equation for a common garden with different nurseries.

We modeled each trait \(y_{bpfr}\) such as:

\[\begin{equation} \begin{aligned} y_{bpfnr} & \sim \mathcal{N}(\mu_{bpfn},\sigma^{2}_{r})\\[3pt] \mu_{bpfn} & = \beta_{0} + B_{b} + P_{p} + F_{f(p)} + N_{n}\\[3pt] \end{aligned} \end{equation}\]

\(\beta_{0}\) is the global intercept. \(B_{b}\) are the block intercepts with a variance \(\sigma^{2}_{B}\). \(P_{p}\) are the population intercepts with a variance \(\sigma^{2}_{P}\). \(F_{f}\) are the family intercepts with a variance \(\sigma^{2}_{F}\). \(N_{n}\) are the nursery intercepts with a variance \(\sigma^{2}_{N}\). \(\sigma^{2}_{r}\) is the residual variance.

The prior of \(\beta_{0}\) was weakly informative and centered around the mean of the observed values for the trait under considered, as follows:

\[\beta_{0} \sim \mathcal{N}(\mu_{y},20)\]

The population, nursery and block intercepts, \(P_{p}\), \(N_{n}\) and \(B_{b}\) were considered normally-distributed with variances \(\sigma^{2}_{P}\), \(\sigma^{2}_{N}\) and \(\sigma^{2}_{B}\), such as:

\[\begin{bmatrix} B_{b} \\ P_{p} \\ N_{n} \end{bmatrix} \sim \mathcal{N}\left(0, \begin{bmatrix}\sigma^{2}_{B}\\[3pt] \sigma^{2}_{P}\\ \sigma^{2}_{N}\\ \end{bmatrix} \right)\\[3pt]\]

The family intercepts \(F_{f(p)}\) were considered to follow some population-specific normal distributions, such as: \[F_{f(p)} \sim \mathcal{N}(0,\sigma^{2}_{F_{p}})\]

where \(\sigma^{2}_{F_{p}}\) are the population-specific variances among families.

To partition the total variance, we parameterize our model so that only the total variance, \(\sigma_{tot}^{2}\) has a prior, such that:

\[\begin{equation} \begin{aligned} \sigma_{tot}^{2} & = \sigma_{r}^{2} + \sigma_{B}^{2} + \overline{\sigma_{F_{p}}}^{2} + \sigma_{P}^{2} + \sigma_{N}^{2}\\[3pt] \sigma_{r} & = \sigma_{tot} \times \sqrt(\pi_{r})\\[3pt] \sigma_{B} & = \sigma_{tot} \times \sqrt(\pi_{B})\\[3pt] \sigma_{P} & = \sigma_{tot} \times \sqrt(\pi_{P})\\[3pt] \sigma_{N} & = \sigma_{tot} \times \sqrt(\pi_{N})\\[3pt] \overline{\sigma_{F_{p}}} & = \sigma_{tot} \times \sqrt(\pi_{F})\\[3pt] \sigma_{tot} & \sim \mathcal{S}^{*}(0,1,3) \end{aligned} \end{equation}\]

where \(\overline{\sigma_{F_{p}}}\) and \(\overline{\sigma_{F_{p}}}^{2}\) are the mean of the population-specific among-family standard deviations (\(\sigma_{F_{p}}\)) and variances (\(\sigma^{2}_{F_{p}}\)), respectively, and \(\sum_{l}^{5}\pi_{l}=1\) (using the simplex function in Stan).

The population-specific among-families standard deviations \(\sigma_{F_{p}}\) follow a log-normal distribution with mean \(\overline{\sigma_{F_{p}}}\) and variance \(\sigma^{2}_{K}\), such as:

\[\begin{equation*} \begin{aligned} \sigma_{F_{p}} & \sim \mathcal{LN}\left(\ln(\overline{\sigma_{F_{p}}})-\frac{\sigma^{2}_{K}}{2},\sigma^{2}_{K}\right)\\[3pt] \sigma_{K} & \sim \exp(1)\\[3pt] \end{aligned} \end{equation*}\]

Stan code of the models:

Code

model_all_pops_different_nurseries <- stan_model(here("scripts/stan_models/StandardQuantGenModel_8.stan"))
model_all_pops_one_nursery <- stan_model(here("scripts/stan_models/StandardQuantGenModel_7.stan"))
model_all_pops_different_nurseries 

S4 class stanmodel 'anon_model' coded as follows:
// Estimating population-specific heritabilities and evolvabilities in the common gardens with trees from different nurseries
data {
  int<lower=1> n;                                                              // Number of observations
  vector[n] y;                                                                 // Response variable
  int<lower=0> n_bloc;                                                         // Number of blocks
  int<lower=0> n_pop;                                                          // Number of populations
  int<lower=0> n_fam;                                                          // Number of families
  int<lower=0> n_nurs;                                                         // Number of nurseries
  int<lower=0, upper=n_bloc> bloc[n];                                          // Blocks
  int<lower=0, upper=n_pop> pop[n];                                            // Populations
  int<lower=0, upper=n_pop> which_pop[n_fam];                                  // Populations
  int<lower=0, upper=n_fam> fam[n];                                            // Families
  int<lower=0, upper=n_nurs> nurs[n];                                          // Nurseries
}
parameters {
  real beta0; // global intercept
  simplex[5] pi;
  real<lower = 0> sigma_tot;
  real<lower = 0> sigma_K;
  vector[n_pop] z_log_sigma_fam;

  vector[n_bloc] z_bloc;
  vector[n_pop] z_pop;
  vector[n_fam] z_fam;
  vector[n_nurs] z_nurs;

}
transformed parameters {
  real R_squared;

  real<lower = 0>  sigma_r;
  real<lower = 0>  sigma_bloc;
  real<lower = 0>  sigma_pop;
  real<lower = 0>  sigma_nurs;

  real mean_sigma_fam;
  vector[n_pop] sigma_fam;

  vector[n_pop] alpha_pop;
  vector[n_fam] alpha_fam;
  vector[n_bloc] alpha_bloc;
  vector[n_nurs] alpha_nurs;

  vector[n] mu; // linear predictor

  // variance partitioning with the simplex pi
  sigma_r = sqrt(pi[1]) * sigma_tot;
  sigma_bloc = sqrt(pi[2]) * sigma_tot;
  sigma_pop = sqrt(pi[3]) * sigma_tot;

  mean_sigma_fam= sqrt(pi[4]) * sigma_tot;
  sigma_fam = exp(log(mean_sigma_fam) - (square(sigma_K)/2)  + z_log_sigma_fam*sigma_K);

  sigma_nurs= sqrt(pi[5]) * sigma_tot;

  alpha_pop = z_pop*sigma_pop;

  for(f in 1:n_fam){
    alpha_fam[f] =  z_fam[f]*sigma_fam[which_pop[f]];
  }

  alpha_bloc = z_bloc*sigma_bloc;
  alpha_nurs = z_nurs*sigma_nurs;

  mu = rep_vector(beta0, n) + alpha_fam[fam] + alpha_bloc[bloc] + alpha_pop[pop] + alpha_nurs[nurs];
  R_squared = 1 - variance(y - mu) / variance(y);
}
model{
  //Priors
  beta0 ~ normal(mean(y),2);
  sigma_tot ~ student_t(3, 0.0, 1.0);

  z_pop ~ std_normal();
  z_bloc ~ std_normal();
  z_fam ~ std_normal();
  z_nurs ~ std_normal();

  z_log_sigma_fam ~ std_normal();
  sigma_K ~ exponential(1);

  // Likelihood
  y ~ normal(mu, sigma_r);
}
generated quantities {
  //Variances
  real<lower=0> sigma2_r;
  real<lower=0> sigma2_pop;
  real<lower=0> sigma2_bloc;
  vector[n_pop] sigma2_fam;
  real<lower=0> sigma2_nurs;

  // Heritability
  vector[n_pop] h2_pop;

  // Evolvability
  vector[n_pop] I_pop;

  // Posterior predictive check
  vector[n] y_rep;

  // Log-Likelihood (for WAIC/loo computations)
  //vector[n] log_lik;

  // Bayesian R2
  real bayes_R2_res; // residual based R2
  real bayes_R2;     // model based R2, i.e.  using modelled (approximate) residual variance

  sigma2_r = square(sigma_r);
  sigma2_pop = square(sigma_pop);
  sigma2_fam = square(sigma_fam);
  sigma2_nurs = square(sigma_nurs);
  for(p in 1:n_pop)  h2_pop[p] = sigma2_fam[p]/(sigma2_r+sigma2_fam[p]);
  for(p in 1:n_pop)  I_pop[p] = sigma2_fam[p]/(mean(y)^2);// should we divide by the mean(y) of each population?
    sigma2_bloc = square(sigma_bloc);
  for(i in 1:n)  {
    y_rep[i] = normal_rng(mu[i], sigma_r);
    //log_lik[i] = normal_lpdf(y[i]| mu[i], sigma_r); // log probability density function
  }

  bayes_R2_res = variance(mu) / (variance(mu) + variance(y-mu));
  bayes_R2 = variance(mu) / (variance(mu) + sigma2_r);

} 

Code

lapply(site_codes, function(site_code){
  
  
  subsite <- subdata %>% 
    filter(site==site_code) %>% 
    arrange(pop)
  
  if(site_code %in% c("FE","FS")){ 
    model <- model_all_pops_one_nursery
    pars <- c("beta0", "pi", 
              "bayes_R2_res", "bayes_R2",
              "h2_pop","I_pop",
              "alpha_bloc","alpha_pop","alpha_fam",
              "sigma2_r","sigma2_bloc","sigma2_pop","sigma2_fam")
  } else if(site_code=="FW"){
    model <- model_all_pops_different_nurseries
    pars <- c("beta0", "pi", 
              "bayes_R2_res", "bayes_R2",
              "h2_pop","I_pop",
              "alpha_bloc","alpha_pop","alpha_fam","alpha_nurs",
              "sigma2_r","sigma2_bloc","sigma2_pop","sigma2_fam","sigma2_nurs")
      
  }
  
  stanlist <- compose_data(subsite)
  
  stanlist$which_pop <- as.numeric(as.factor(pull(unique(subsite[c("pop","fam")])[,"pop"])))
  
  stanfit <- sampling(model,
                      data = stanlist,
                      pars=pars,
                      iter = n_iter, 
                      chains = n_chains, 
                      cores = n_chains,
                      save_warmup = save_warmup,
                      thin=n_thin)
}) %>% 
  setNames(site_codes) %>% 
  saveRDS(file=here(paste0("outputs/Evolvability/",trait,"_FullModel.rds")))

Warning: There were 68 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: There were 72 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Warning: There were 8 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

The models have divergent transitions. We can visualize them with pairs plots.

Code

readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel.rds"))) %>% 
lapply(function(x){
 mcmc_pairs(as.array(x), 
           np = nuts_params(x), 
           pars = c("pi[1]","pi[2]","pi[3]","pi[4]"),
           off_diag_args = list(size = 0.75)) 
  
})

$FE

$FW

$FS

Comparing \(h^2\) and \(I\) estimates

We check the correlations among \(h^2\) and \(I\) values estimated calculated with population-specific models or with the full model including all populations.

Code

list_df <- lapply(c("I","h2"), function(x){
  

# Extract the median estimates of h2/I in the population-specific models
# ======================================================================
df_pop_specific_models <- readRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels.rds"))) %>% 
  lapply(function(site_i){
    
    site_i %>% lapply(function(pop_i){
    
  broom.mixed::tidyMCMC(pop_i,pars=(x),
                droppars = NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) 
      
    }) %>% bind_rows(.id = "pop")
  }) %>% bind_rows(.id = "site") %>% 
  mutate(model="pop_specific_models") %>% 
  mutate(model_site=paste0(model,"_",site)) %>% 
  dplyr::select(-term)

# Extract the median estimates of h2/I in the population-specific models
# ======================================================================
df <- readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel.rds"))) %>%
  lapply(function(site_i){
    
  broom.mixed::tidyMCMC(site_i,pars=(paste0(x,"_pop")),
                droppars = NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) 
      
  }) %>% 
  bind_rows(.id = "site") %>% 
  mutate(model="full_model") %>% 
  mutate(model_site=paste0(model,"_",site)) %>% 
  mutate(pop=rep(unique(subdata$pop) %>% sort,3)) %>% 
  dplyr::select(-term) %>% 
  bind_rows(df_pop_specific_models)

  
}) %>% setNames(c("I","h2"))


list_df[["h2"]] %>%   
  ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site,shape=model)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  xlab("") +
  ylab("Narrow-sense heritability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank())  +
    guides(color=guide_legend(ncol=1))

Code

list_df[["I"]] %>%   
  ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site,shape=model)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  xlab("") +
  ylab("Evolvability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank())  +
    guides(color=guide_legend(ncol=1))

Code

list_df[["I"]] %>% 
  dplyr::select(estimate,model_site,pop) %>% 
  pivot_wider(names_from=model_site,values_from=estimate) %>% 
  dplyr::select(-pop) %>% 
  cor() %>% 
  corrplot(method = 'number',
           type = 'lower', diag = FALSE,
           mar=c(0,0,2,0),
           title="Evolvability",
           number.cex=0.9,tl.cex=0.9)

Code

list_df[["h2"]] %>% 
  dplyr::select(estimate,model_site,pop) %>% 
  pivot_wider(names_from=model_site,values_from=estimate) %>% 
  dplyr::select(-pop) %>% 
  cor() %>% 
  corrplot(method = 'number',
           type = 'lower', diag = FALSE,
           mar=c(0,0,2,0),
           title="Narrow-sense heritability",
           number.cex=0.9,tl.cex=0.9)

Population-specific models - filtered datasets

One potential explanation for the divergent transitions is the insufficient number of individuals in some families or families in populations.

Sample size

We look at the number of individuals in each family and we count the number of families with 1, 2, 3 or 4 individuals:

Code

nb_ind_per_family <- lapply(site_codes, function(site_code){
  
  subdata %>%
    filter(site == site_code) %>% 
    group_by(fam) %>% 
    summarise(nb_individual_per_family=n())  %>% 
    group_by(nb_individual_per_family) %>% 
    summarise(nb_family=n()) %>% 
    setNames(c("Number of individuals per family", "Number of families"))
}) %>% setNames(site_codes)

In FE:

Code

nb_ind_per_family[["FE"]] %>%  kable_mydf

Number of individuals per family	Number of families
1	1
2	3
3	13
4	152

In FS:

Code

nb_ind_per_family[["FS"]] %>%  kable_mydf

Number of individuals per family	Number of families
2	2
3	15
4	151

We also look at the number of families per population:

Code

nb_fam_per_pop <- lapply(site_codes, function(site_code){
  
  subdata %>%
    filter(site == site_code) %>% 
    dplyr::select(pop,fam) %>% 
    distinct() %>% 
    group_by(pop) %>% 
    summarise(nb_fam_per_pop=n()) %>% 
    setNames(c("Population", "Number of families per population"))
}) %>% setNames(site_codes)

In FE:

Code

nb_fam_per_pop[["FE"]] %>%  kable_mydf

Population	Number of families per population
AB	8
AC	8
AM	8
BB	8
BE	9
BW	8
CC	8
CG	8
CR	8
GA	8
GC	8
GD	8
GE	8
GL	8
GT	8
LC	8
MG	8
RD	8
RM	8
SD	8
SO	8

In FS:

Code

nb_fam_per_pop[["FS"]] %>%  kable_mydf

Population	Number of families per population
AB	8
AC	8
AM	8
BB	8
BE	8
BW	8
CC	8
CG	8
CR	8
GA	8
GC	8
GD	8
GE	8
GL	8
GT	8
LC	8
MG	8
RD	8
RM	8
SD	8
SO	8

Keeping only families with 4 individuals

We fit again the models with only families in which there are 4 individuals, so only in FE and FS.

We identify the families with four individuals:

Code

site_codes <- c("FE","FS")

selected_families <- lapply(site_codes, function(site_code){
  
  subdata %>%
    filter(site == site_code) %>% 
    group_by(fam) %>% 
    summarise(nb_individual_per_family=n()) %>% 
    filter(nb_individual_per_family==4) %>% 
    pull(fam)
}) %>% setNames(site_codes)

We look at the trait distribution after keeping only populations with 4 individuals.

Code

  subdata %>%
    filter(site %in% site_codes) %>% 
    filter((site == "FE" & fam %in% selected_families[["FE"]]) |
           (site == "FS" & fam %in% selected_families[["FS"]])) %>% 
    ggplot(aes(x=y)) +  
    geom_histogram(aes(y=after_stat(density)), colour="blue",fill="white",bins = 34) +
    geom_density(alpha=.2,fill="pink") +
    ggtitle(paste0(trait," distribution after keeping only populations with 4 individuals")) +
    facet_wrap(~site,scales="free") + 
    theme_bw() 

We run the population-specific models and the full models:

Code

lapply(site_codes, function(site_code){
  
  subdata %>%
    filter(site == site_code) %>% 
    filter(fam %in% selected_families[[site_code]]) %>% 
    dplyr::select(pop,fam) %>% 
    distinct() %>% 
    group_by(pop) %>% 
    summarise(nb_fam_per_pop=n()) %>% 
    group_by(nb_fam_per_pop) %>% 
    summarise(nb_family=n())
}) %>% setNames(site_codes)

$FE
# A tibble: 4 × 2
  nb_fam_per_pop nb_family
           <int>     <int>
1              5         1
2              6         4
3              7         5
4              8        11

$FS
# A tibble: 4 × 2
  nb_fam_per_pop nb_family
           <int>     <int>
1              5         1
2              6         3
3              7         8
4              8         9

Code

  # Population-specific models
  # ==========================
lapply(site_codes, function(site_code){
  
 subsite <- subdata %>% 
    filter(site==site_code) %>% 
    filter(fam %in% selected_families[[site_code]]) %>% 
    dplyr::select(bloc, pop, fam, y) %>% 
    arrange(pop)
  
 lapply(unique(subsite$pop), function(pop_i){
   
   subpop <- subsite %>% 
    filter(pop==pop_i) 
   
    model <- model_pop_specific_one_nursery
    pars <- c("beta0", "pi", "R_squared", "h2","I",
              "alpha_block","alpha_fam",
              "sigma2_r","sigma2_block","sigma2_fam","sigma2_tot")
 
  stanfit <- sampling(model,
                      data = compose_data(subpop),
                      pars=pars,
                      iter = n_iter, 
                      chains = n_chains, 
                      cores = n_chains,
                      save_warmup = save_warmup,
                      thin=n_thin)   
 }) %>% 
   setNames(unique(subsite$pop))
}) %>% 
    setNames(site_codes) %>% 
    saveRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels_FilteredDataset.rds")))
  
  
    # Full model
    # ==========
  
lapply(site_codes, function(site_code){
  
  
  subsite <- subdata %>% 
    filter(site==site_code) %>% 
    filter(fam %in% selected_families[[site_code]]) %>% 
    dplyr::select(bloc, pop, fam, y) %>% 
    arrange(pop)
  
    model <- model_all_pops_one_nursery
    pars <- c("beta0", "pi", 
              "bayes_R2_res", "bayes_R2",
              "h2_pop","I_pop",
              "alpha_bloc","alpha_pop","alpha_fam",
              "sigma2_r","sigma2_bloc","sigma2_pop","sigma2_fam")

  
  stanlist <- compose_data(subsite)
  
  stanlist$which_pop <- as.numeric(as.factor(pull(unique(subsite[c("pop","fam")])[,"pop"])))
  
  stanfit <- sampling(model,
                      data = stanlist,
                      pars=pars,
                      iter = n_iter, 
                      chains = n_chains, 
                      cores = n_chains,
                      save_warmup = save_warmup,
                      thin=n_thin)
}) %>% 
  setNames(site_codes) %>% 
  saveRDS(file=here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds")))

Warning: There were 99 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: Examine the pairs() plot to diagnose sampling problems

Code

readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds"))) %>% 
lapply(function(x){
 mcmc_pairs(as.array(x), 
           np = nuts_params(x), 
           pars = c("pi[1]","pi[2]","pi[3]","pi[4]"),
           off_diag_args = list(size = 0.75)) 
  
})

$FE

$FS

Code

readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds"))) %>% 
lapply(function(x){
 mcmc_pairs(as.array(x), 
           np = nuts_params(x), 
           pars = c("alpha_pop[2]","h2_pop[2]","h2_pop[17]","sigma2_r"),
           off_diag_args = list(size = 0.75)) 
  
})

$FE

$FS

Comparing estimates

We check the correlations between \(h^2\) and \(I\) estimated population-specific models or with the full model including all populations.

Code

list_df <- lapply(c("I","h2"), function(x){
  

# Extract the median estimates of h2/I in the population-specific models
# ======================================================================
df_pop_specific_models <- readRDS(file=here(paste0("outputs/Evolvability/",trait,"_PopSpecificModels_FilteredDataset.rds"))) %>% 
  lapply(function(site_i){
    
    site_i %>% lapply(function(pop_i){
    
  broom.mixed::tidyMCMC(pop_i,pars=(x),
                droppars = NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) 
      
    }) %>% bind_rows(.id = "pop")
  }) %>% bind_rows(.id = "site") %>% 
  mutate(model="pop_specific_models") %>% 
  mutate(model_site=paste0(model,"_",site)) %>% 
  dplyr::select(-term)

# Extract the median estimates of h2/I in the population-specific models
# ======================================================================
df <- readRDS(here(paste0("outputs/Evolvability/",trait,"_FullModel_FilteredDataset.rds"))) %>%
  lapply(function(site_i){
    
  broom.mixed::tidyMCMC(site_i,pars=(paste0(x,"_pop")),
                droppars = NULL, ess = F, rhat = F, conf.int = T,conf.level = conf_level) 
      
  }) %>% 
  bind_rows(.id = "site") %>% 
  mutate(model="full_model") %>% 
  mutate(model_site=paste0(model,"_",site)) %>% 
  mutate(pop=rep(unique(subdata$pop) %>% sort,2)) %>% 
  dplyr::select(-term) %>% 
  bind_rows(df_pop_specific_models)

  
}) %>% setNames(c("I","h2"))


list_df[["h2"]] %>%   
  ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site,shape=model)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  xlab("") +
  ylab("Narrow-sense heritability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank())  +
    guides(color=guide_legend(ncol=1))

Code

list_df[["I"]] %>%   
  ggplot(aes(y = estimate, x = pop, ymin = conf.low, ymax = conf.high,color=site,shape=model)) +
  geom_pointinterval(position = position_dodge(width = .4),point_size=3,alpha=0.6,size=5) +
  xlab("") +
  ylab("Evolvability") +
  theme(axis.text = element_text(size=18),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank())  +
    guides(color=guide_legend(ncol=1))

Code

list_df[["I"]] %>% 
  dplyr::select(estimate,model_site,pop) %>% 
  pivot_wider(names_from=model_site,values_from=estimate) %>% 
  dplyr::select(-pop) %>% 
  cor() %>% 
  corrplot(method = 'number',
           type = 'lower', diag = FALSE,
           mar=c(0,0,2,0),
           title="Evolvability",
           number.cex=0.9,tl.cex=0.9)

Code

list_df[["h2"]] %>% 
  dplyr::select(estimate,model_site,pop) %>% 
  pivot_wider(names_from=model_site,values_from=estimate) %>% 
  dplyr::select(-pop) %>% 
  cor() %>% 
  corrplot(method = 'number',
           type = 'lower', diag = FALSE,
           mar=c(0,0,2,0),
           title="Narrow-sense heritability",
           number.cex=0.9,tl.cex=0.9)

References

Archambeau, Juliette, Marta Benito Garzón, Marina de Miguel, Benjamin Brachi, Frédéric Barraquand, and Santiago C González-Martı́nez. 2023. “Reduced Within-Population Quantitative Genetic Variation Is Associated with Climate Harshness in Maritime Pine.” Heredity, 1–11.

Beaton, Joan, Annika Perry, Joan Cottrell, Glenn Iason, Jenni Stockan, and Stephen Cavers. 2022. “Phenotypic Trait Variation in a Long-Term Multisite Common Garden Experiment of Scots Pine in Scotland.” Scientific Data 9 (1): 671.